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Context. The current star formation models imply that the binary fraction of population III stars is non zero. The evolution of such 
binaries must have led to formation of compact object binaries. 

Aims. In this paper we estimate the gravitational wave background originating in such binaries and discuss its observability . 
Methods. The properties of the population III binaries are investigated using a binary population synthesis code. We numerically 
model the background and we take into account the evolution of eccentric binaries. 

Results. The gravitational wave background from population III binaries dominates the spectrum below 100 Hz. If the binary fraction 
is larger than 10~ 2 the background will be detectable by LISA and DECIGO. 

Conclusions. Gravitational wave background from population III binaries will dominate the spectrum below 100 Hz. LISA, ET and 
DECIGO should either see it easily or, in case of non detection, provide very strong constraints on the properties of the population III 
stars. 
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Coalescing compact object binaries are the most promising 
sources of gravitational waves. This include the double neutron 
star binaries (DNS), the mixed, black hole neutron star bina- 
ries (BHNS) and binary black holes (BBH). Gravitational wave 
signal may be detected as single events from nearby sources, 
however for the distant ones the signals will overlap contribut- 
ing to gravitational wave background (GWB). GWB from coa- 
lescing compact obj ect binaries has been investigated by numer- 
ous authors (e . g..lTu mhnson & Shulfll2000t[Bromm et al.ll200lt 
Rosadd I201 U IZhu et alj |201 U lOlmez et all |201 lflRegimbau 



201 it iMarassi et al.ll2.QlTh . The existing gravitational wave de- 



tectors LIGO and VIRGO have achieved their design sensi- 
tivity (lAccadia et al.1 d2011l) . lKawabe & the LIGO Collaboration! 
(2008)). Currently they are undergoing transitions to the 
Advanced phase with the sensitivity increased by a factor 30. 
Additionally there a re plans to construct other large scale d etec- 
tors such as LCGT dKuroda & LCGT Collaborationll20irj) . In a 
more distant fut ure third generation d etector like the Einstein 
Telescope (ET) dVan Den BroecklfeOlO ) should be constructed. 
Theses instrumental developments make the investigation of 
stochastic backgrounds of gravitation waves very important. We 
expect that at lower frequencies such as designed band of ET (~ 
1Hz), large population of distant object will be observed. The 
biggest contribution in that frequency window will come from 
massive BBH. There are no direct observations of that kind of 
objects, therefore there are many uncertainties in stellar evo- 
lution models leading to creation BBH. Because of the range 
of future detectors there is possibility that among observed ob- 
jects could be the remnants of the oldest stars in the Universe. 
Population III stars should be metal free, therefore they should 
produce black holes of masses reaching hundreds solar masses. 
If our understanding of their evolution is correct, there should 



be significant gravitational background from remnants of popu- 
lation III stars. 



In this paper we investigate the possible contribution to the 
GWB from population III star binaries. Population III stars are 
the first stars in the Universe that formed out of metal free 
material. As such, they had very different properties than the 
currently forming stars. The evolution and prop erties of sin- 
gle Po p ulation III star have been investigated by Baraffe et al. 
(120011) . QUI d2008l) . They found that the single population III 
star are stable up to a few hundred solar masses. Moreover they 
evolve with a little mass loss. The evolution of such massive 
star will end up either in a pair instability supernova leaving 
no remnant, or forming a black hole through a direct collapse. 
While no population III stars are known, their properties have 
been investigated by numerical simulations of the collapse of 
metal free clouds of gas. These simulations have shown that 
the initial mass function of population III star is skewed to- 
wards higher masses than the one on the population I or II stars 
(lAbel et alJl200llBrommer al. 2002). The existence of binaries 
among Population III star has been uncertain, however simula- 
tions of collapse of rotating clouds of gas h ave shown bar in- 
stabilities and form ation of two protostars (Saigo et al. 2004; 
iMachida eta"DT 2008). Since all known stellar populations con- 
tain a significant fraction of binaries, it seems justified to con- 
sider the properties of Population III binaries. The evolution of 
Popula tion III binaries has been in vestigated by Belczvnski et al. 
(2004) and Kulczycki et al. (2006). They found that coalescence 
of BBH binaries originating in population III stars may be a sig- 
nificant source of gravitational waves for the current and future 
gravitational wave interferometers. In section 2 we analyze the 
gravitational wave spectrum form eccentric binaries, section 3 is 
devoted to the analysis of properties of population III binaries. 
Section 4 contains an outline of the estimate of the GWB. We 
present the results in section 5 and in section 6 we discuss them. 
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2. Population of Pop. Ill stars - models 

The evolution of population III metal free stars was firs t 
examined in detail about ten years ag o iHeger et all (|2001 | ) . 
Bara ffeet all (120011) . iMarigo et all (l200lb . iMarigoet all (|2003), 
Schaerer (200^Tlt was found that such stars are stable with 
masses up to 500 M . Population III star evolve very quickly 
with nearly no mass loss. Depending on the initial mass they 
may form neutrons stars, black holes or leave no remnant at 
all, because of total star disruption in pair instability supernova 
(Hege r&Wooslevl l2002). The mass spectrum of population III 
stars is still a matter of a debate. On one hand the simulations 
of metal free cloud collapse indicate the possibility of forming 
high mass objects. The searches for low mass metal free stars 
that should still be around brought no success, which indicates 
that the low mass cutoff of population III initial mass function 
was significantly higher than for the current population of stars. 
An additional question is whether population III stars also were 
formed as binaries. Every known stellar population contains a 
significant fraction of binaries, so one should also expect this to 
be the case for population III stars. Numerical simulations of col- 
lapse of rotating metal free clouds indicate appearance of a bar 
instability which le ads to formation of two pro tostars in a binary 
( Saig o et al.1 120041). In a recent simulations dTurket al.1 [2 009; 



Initial parameters of compact binaries 



Stacv et al.l2010HGreif et alJ201 UlClark et al.l201 lHStacv et alJ 
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2012) a similar result was shown. A massive protostellar cloud 



develops spiral structure, instabilities and leads to formation of 
several multiple systems. Thus we conclude that there is ample 
evidence to allow investigation of population III binary evolution 
and it consequences. 

We model the evoluti on of population III binaries us- 
ing a code described in iBelczvnski et alJ d2004l). The cod e 
uses the evolutionary sequences b a sed on IHeger et ail d2001l) . 
iBaraffe et all d2001l) , iMarigo et all d2001l) . It includes a detail 
treatment of the calculation of stellar remnants. The binary evo- 
lution calculation involves analysis of stability of mass trans- 
fers. The common envelope mass transfers end in a merger if 
the donor is a main se quence star and for giant stars we use the 
standard formalism o f lWebbinkl(fT98l . 

We assume that the initial mass function of the primaries 
stretches between 10 and 500 M with the exponent of -2.35. 
The mass ratio is drawn from a flat distribution. The initial or- 
bital separation distribution is flat in logarithm of a. The range of 
initial orbital separation is bounded from below by the require- 
ment that the stars are not in contact at zero age main sequence 
(ZAMS) - a m j„ = 13(R\ + Rz)- The upper bo und is more ar- 
bitrary . We have chosen a max = \Q 6 R & . While Bro mm & LoebJ 
(2006) argue that the upper limit on initial separation is as large 
as 5 x 10 7 ^ o , we wanted to include all systems that interact. The 
binaries with initial separations close to or above our upper limit 
will not interact and if they survive the evolution and form com- 
pact object binaries, they will only contribute to the GW back- 
ground spectrum below 10~ n Hz. The uncertainty in the number 
of binaries due to variation of the upper cutoff a max is small be- 
cause of logarithmic distribution of initial orbital separations. 
We calculate the evolution of binaries and note that cases when 
a compact object binary was formed. Then we trace it evolution 
due to emission of gravitational waves. 



2. 1 . Properties of compact objects 

We evolved 10 6 binary stars of zero metallicity using popula- 
tion synthesis code Star/Track. As a result we obtained 462496 
compact objects, which we treat as a realistic sample of pop- 
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Fig. 1. The initial properties of the Population III compact ob- 
ject binaries. The top panel shows the distribution of the orbital 
properties, ie. the density in the plane spanned by the semi-major 
axis and the eccentricity of the orbit. The lower panel shows the 
distribution of the chirp masses and the merger times of the bi- 
naries. Vertical line represents Hubble time. Numbers represent 
percentage of all binaries in simulation. 

ulation III remnants. More than 60% of them are BBH sys- 
tems, where both components have masses greater than 2.5M Q . 
Detailed distributions of initial parameters are presented on 
Figure I2TT1 On the top panel one can see that there is very broad 
range of initial eccentricities - over 60% of all binaries have ec- 
centricity greater than 0.1. Bottom panel shows distributions of 
chirp masses and coalescence times. Only 12% of binaries will 
merge within present Hubble time. 

3. Gravitational wave spectrum of eccentric 
binaries 

While the spectrum from circular binaries is well known, we 
have shown above that a large fraction of population III binaries 
will be eccentric. It is therefore important to investigate the pos- 
sible effects of eccentrici ty on the gravitational wave spectra o f 
compact object binaries. (Peters & Mathews 1 963t iPetersi ri 964 ) 

3.1. General derivation 

The gravitational wave spectrum from a binary can be calculated 

as: 
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dfgw 
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(1) 



The velocity on an eccentric orbit is changing over its period, 
therefore the instantaneous orbital frequency is also varying. 
Thus the system is radiating in some range of frequencies, not 
just in one particular, like in the case of a circular orbit. 

In the quadrupole approximation the orbit decays due to 
gravitational wave emission and the orbital parameters change 
as: 
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while orbital frequency is 

forb 



1 /GM\" 2 



(4) 



where we denoted the total mass of the binary M, jj is the re- 
duced mass, and a is the semi-major axis. An eccentric binary 
emits gravitational waves in a spectrum of harmonics of the or- 
bital frequency 



few ~ n forb 7 



dfgn 
dt 



df 



orb 



dt 



(5) 
(6) 



In the case of circular orbit gravitational waves are emitted only 
in one mode n — 2. Inserting equations [4] and [2] into equation [6] 
we obtain: 



df% 96 (2x\ >13 r^G^M* 

— c nc), 

n ) c 3 



dt 



(7) 



where we introduced the chirp mass, which is quantity that deter- 
mines amplitude and frequency dependence of the inspiral grav- 
itational wave sig nal. M chi rp = yU 3/s M 2/s . The power of radiatio n 
for each harmonic was calculated bv lPeters & Mathews! (1 19631) : 



dE , 32 G 4 fi 2 M 3 , 

* (n) = T^5-* (n ' e) ' 



(8) 



where g(n, e) is: 
g{n,e) 



rr I 2 
~ \ Un-iine) - 2eJ n -i(ne) + -J„(ne) 
32 { n 

+2eJ n+ \{ne) - J„+ 2 (ne)] 2 

+(1 - e 2 ) [J„- 2 (ne) - 2J„{ne) + J„+ 2 {ne)] 2 



in 1 



(9) 



and J„ are the Bessel functio ns. This result follows fro m Fourier 
analysis of the Kepler motion [Peters & Mathews! (1 19631) . We ob- 
tain the instantaneous spectrum of gravitational waves from an 
eccentric binary by combining equations [BE] an d 13 
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In the case of the circular orbit this reduces to the well known 
result, e.g. (lPhinnevl2001l) : 
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since ^(e = 0) = 1, g(n, e — 0) = 1 when n — 2 and g(n + 2, e — 
0) = 0. The result has to be averaged over the lifetime of the 
binary. In the following section we discuss two cases: the long 
lived and the short lived binaries, as outlined below. 

3.2. Long-lived systems 

By term long-living we mean that the time to coalescence is 
longer than present Hubble time, which we assume for simplicity 
tn = 10 4 Myr. We are interested in objects that have not merged. 
Such objects have very wide orbits and the semi-major axis of 




Fig. 2. Spectra of the GW radiation from long-living binary 
(Mi = M 2 = 300 M a = 1.2 x 10 6 fl o ) for different values of 
eccentricity. A single 'plus' mark corresponds to circular orbit, 
by crosses we show case for eccentricity 0.5, stars corresponds 
to eccentricity 0.7 and squares to initial eccentricity 0.9 



the binary as well as its eccentricity can be assumed to be con- 
stant during the entire evolution time. 

In this case the spectrum is discrete as the binary emits the 
gravitational radiation only at the frequencies specified by the 
harmonics. 



dE V xr f t \ dE 
Tr ^f-nforb)— ; 



(12) 



where the spectrum in each harmonic is given by equationfTOl 

To illustrate this case we chose binary consisting of two 
massive black holes: M\ = M 2 = 300 M on wide orbit a - 
1.2 x 1O 6 ^ . We present the spectra for a few values of the ec- 
centricity in Figure|2] 

When eccentricity is almost zero we have just one point - it 
corresponds to very wide, circular orbit. Velocity in that case is 
constant and we observe only one frequency. While eccentricity 
is increasing we observe that spectrum is wider. Also maximum 
of energy distribution is shifted to the higher frequencies. For the 
highest considered eccentricity (0.9), there is a minimum around 
frequency 10" 11 . It is caused by the specific structure of the func- 
tion described by equation|9] 

3.3. Short-lived systems 

In that group we consider systems which can coalescence dur- 
ing they evolution. The orbit is changing and we observe ra- 
diation in very wide range of frequencies. In this case the bi- 
nary emits a continues spectrum over its lifetime. The orbital 
frequency changes from //,„ - the initial orbital frequency deter- 
mined by the orbital parameters after formation of the second 
compact object in the binary, to the final orbital frequency just 
before merger. For the final orbital frequency we chose the fre- 
quency at the marginally stable orbit 
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As the binary evolves, its eccentricity decreases and the evolu- 
tion of eccentricity can be found by solving the equation: 
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Fig. 3. Spectrum of the GW radiation from short-lived binaries 
for different values of the initial eccentricity. Solid line corre- 
sponds to circular orbit, by dashed line we show case for eccen- 
tricity 0.5, dotted line corresponds to eccentricity 0.7 and dash- 
dotted line to initial eccentricity 0.9 



Let us denote the dependence of the eccentricity on the orbital 
frequency as e(f or b). 

The gravitational wave spectrum is 



dE dE 
df Z-i ,/ /<; "" 



dfg W 



where the functions 4E- are evaluated only for the frequencies 

in the range nf ini < f' g l w < nf lsco . 

We illustrate this case with a binary with the same mass as in 
previous section: M\ - = 300 M , but a much tighter orbit 
a - 1 .2 x 10 2 R Q . We present the spectra for four different values 
of eccentricity in Figure [3] 

For the case of zero eccentricity the spectrum is a simple 
power law. For non zero eccentricities the power law spectrum 
acquires a tail that increases to higher frequencies. The shape of 
the high frequency tail in the spectrum is determined by value of 
the eccentricity of the system. 

4. Gravitational wave background radiation 

We analyze population of the massive black hole binaries us- 
ing th e Star/Track population synthesis code (Belczvnski et al.l 
|2002|) . For each system we obtain its component masses Mij, 
M%,i, and its initial orbital parameters: the semi-major axis a,-, 
and the eccentricity e,. We denote the total mass of the system 



as Mi - M\j + M2J , the reduced mass as /j, 
chirp mass M chirpJ = jt? /s M? /5 . 



M Li M 2j 



and the 



4.1. Evolution of the orbit 

Initial value of semi-major axis is given by synthesis population 
code and we can use it to estimate the lowest frequency possible 
for particular binary: 



Ai = 



GM 



^n 2 



(16) 



Eccentricity and semi-major axis will be changing during the 
evolution according to two differential equations: 
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(18) 



We have to calculate the size of the orbit after Hubble time. 
In case of short-living systems we assume final great semi-major 
axis as aj — 10~ 3 7? o (beyond that value we should consider tidal 
effects). Angular momentum time scale is comparable with age 
of the binary, so it will merge eventually. In the other case (long- 
living systems) we have to solve equations ( TTTb and (fT8l . Then, 
we can calculate frequency corresponding to the final size of the 
orbit in the same way as we did it for initial parameters. 

4.2. Spectrum 

A stochastic background is described in terms of present-day 
gravitational wave energy density (p gw ) per logarithmic fre- 
quency int erval normaliz ed to the critical rest-mass energy den- 
sity (p,.c 2 ) lPhinnevl(l2001l) : 



(15) n gw (f) = 



1 dp gw (f) 
p c c 2 d\n(f) 



(19) 



The total present day energy density contained in the gravita- 
tional radiation is related to Cl(f) and amplitude of the grav- 
itational wave spectrum over a logarithmic frequency interval 
(h c {f)) through: 

f 00 9 df f °° nc 2 , , df 

£ gw = J o p c c 2 n gw (f)j= J o 4^/^(/)j- (20) 

On the other hand & gw should be the sum of the energy radiated 
by all sources at all redshifts. 
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Since Universe is expanding, the observed frequency is different 
from the emitted and it is equal to f r = /(l + z), while N(z) is 
the number of sources in the interval (z,z+ 1). 
Combining ( f20b and (f2Tb we obtain: 
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Using (fTUb we obtain equations describing spectrum from a sin- 
gle harmonic for one type of source: 
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where iVo is the density of one type of sources: 
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N(z) 



-dz. 



The term in the brackets is: 

< <1+d ~" J >4£"a + ;>« 

In the case of Population III binaries the right hand side inte- 
gral is trivial, since we assumed that all Population III stars were 
formed at single redshift z = 15. Taking into account very short 
evolution timescales of massive stars Population III stars, we can 
assume that the binary compact objects were formed at the same 
redshift. 

In order to obtain the energy density from all simulated ob- 
jects, we sum over harmonics (index n) and over the total number 
of binaries (index ;'): 



N 

i=l n=2 3n PcC 

x(/f) 2/3 ^Vo((l+z)- 1/3 ). 
x r(e) x ' 

The amplitude of the gravitational wave spectrum (h c (f)) and 
the spectral density of the gravitational wave background (S (/)) 
follow: 



clnrpj 



Table 1. List of models used in the parameter study 

(25) Model Description 

A Standard 

Bl Binary fraction decreased to 1CT 2 

B2 Binary fraction decreased to 10" 3 

C Star masses drawn independently from same distribution 

D 1 The minimum mass increased to 50M o 

D2 The maximum mass deceased to 100M o 

El The IMF exponent decreased to 1.5 

E2 The IMF exponent increased to 3 

Fl Formation redshift Zf orm = 10 

F2 Formation redshift Zform = 20 

G All binaries are initially circular 



(26) 



h 2 Af) 



^W) 

nf 2 



= 7' 



(27) 



5. Results 



We present the resulting gravitational wave background spec- 
trum in Figure |U The standard model spectrum is shown as the 
thick continuous line. The shape of the background spectrum 
is determined by three different factors. In the low frequency 
regime, below 10~ 5 Hz , the background originates in long lived 
systems, so the spectrum is basically determined by the distribu- 
tion of initial parameters of the binaries at the time of formation. 
In the intermediate region, between 10 4 Hz and 50Hz the spec- 
trum is due to the short lived systems that merger within the 
Hubble time. It has a typical slope of I, which corresponds to 
the orbit decay due to gravitational wave emission. At the fre- 
quencies above 50Hz the spectrum is determined by the higher 
harmonics arising due to eccentricity of the binaries. Note that 
in this work we neglect the merger and ringdown phases which 
may alter the shape of the background spectrum in the region 
above 50Hz. 



5. 1 . Dependence on the models 

The calculation of the background involved assuming the values 
of several parameters. In this section we present the dependence 
of the final result on the choice of those parameters. The list of 
models is shown in TableQ] In the models Bl, B2, C, Dl and D2 
we vary the initial mass distributions of the binaries. In model Dl 
the we have increased the minimum mass of a population III star, 
while in model D2 we decreased upper end of IMF. Models Bl 
and B2 have different values of the initial mass function slope, 
and in model C we draw the two masses in the binary form the 
same IMF. Models El and E2 are characterized by a decreased 
value of the binary fraction to 0.01 for the first one and 0.001 
for the latter one. In models Fl and F2 we place the binaries at 
different formation redshifts, Zform = 10 and 20. Finally model 
G corresponds to the population of circular binaries. 

The shape of the gravitational wave background calculated 
with these different assumptions is shown in Figure [4] The ef- 
fects if decreasing the binary fraction are simply that of the nor- 
malization of the spectrum without altering its shape. Varying 



the initial mass function can change both the normalization of 
the spectrum, as well as the shape of the high frequency region. 
Model D2 has cutoff at frequency 10~ 6 Hz due to the lack of high 
mass binaries. Drawing the two stars from the same distribution 
leads to formation of a large number of small mass binaries and 
the overall level of the spectrum decreases. This is similar to in- 
creasing the IMF exponent, see model E2. In the case when the 
IMF exponent is increased (model El) the level of the spectrum 
goes up as there are more high mass binaries. The gravitational 
wave spectra depend very weakly on the redshift of formation of 
population III stars, as is clearly seen from equation [27] Finally 
the effects of eccentricity are demonstrated by model G where all 
binaries all initially circular. This affects only the high frequency 
tail of the spectrum. 

It is also intersting to compare the gravitational wave back- 
ground from population III binaries with the gravitational wave 
backgrounds from other types of comact binaries. Recently 
Rosado| (1201 ll) performed such a comprehensive calculation. We 
present t he sta ndard model results along with the results of the 
iRosadol d201 ll) calculation in Figure [5] The lower mass comact 
object binaries dominate the region above 100Hz. In the region 
below 100Hz the spectrum of the population III background is 
an order o f magnitude hig her than the NSNS background as cal- 
culated bv lRosadol (1201 lh . 

From Figure|4]one can see that the gravitational wave back- 
ground should be clearly detectable with the next generation in- 
struments. ET will detect it, if the binary fraction was larger than 
10~ 2 . In the case of LISA, the background would show up just 
above the threshold in the band around 10~ 3 Hz. However, for 
future experiments like DECIGO, the gravitational wave back- 
ground from population III stars could be the dominating noise 
source even if the binary was as low as 10 3 . However, at this 
level the gravitational wave background from population III stars 
will be well buried under the background coming from double 
neutron star binaries formed at later epochs. 



6. Summary 

We have analyzed the population III metal free binaries. The 
evolution of such binaries lead to formation of binary black 
holes, which in turn will be a source of gravitational waves. We 
calculate the stochastic background from such binaries and show 
that it can potentially be a significant contribution to the overall 
gravitational wave background. The spectrum has a characteris- 
tic slope of | below a 50Hz, and is falling of above that fre- 
quency. 

We have investigated dependence of the shape of the grav- 
itational wave spectrum on the evolutionary parameters of the 
population III binaries. The evolutionary parameters mainly af- 
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feet the level of the spectrum, while its shape is rather insensitive Van Den Broeck, C. 2010, ArXiv e-prints 
to the parametrization. Only in the region above 50Hz, where Webbink, R. F. 1984, ApJ, 277, 355 

the spectrum fall off rapidly, the shape of the spectrum varies zhu ' X " J - Howe11 ' E - Re § imbau > T -> m ™, D., & Zhu, Z.-H. 201 1, A P J, 739, 86 
with the model. In that spectral range the shape is mainly af- 
fected by the initial mass function of the population III binaries. 
For the binary fraction greater than 10~ 2 the gravitational wave 
background below 100Hz is dominated by the contribution from 
population III BBH binaries. 

The stochastic background from population III BBH bina- 
ries should be detectable by the next generation instruments like 
LISA, DECIGO, and ET. However, the distinctive feature of this 
background is the break in the region of 50 - 100Hz. This is 
unfortunately below the sensitivity of ET. Thus, it will be very 
difficult to distinguish the origin of the background if detected 
by any of the above mentioned instruments. A detailed inves- 
tigation of the spectrum of the gravitational wave background 
above 50Hz would be most interesting. It could potentially re- 
veal the nature of the sources of gravitational wave background, 
their mass spectrum, as well as show some effects connected 
with eccentricity of the binaries. On the other hand lack of detec- 
tion of the stochastic background will lead to some constraints 
on the binary fraction and masses of population III binaries. 
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Fig. 4. The gravitational wave background form population III stars. The standard model results are presented as the thick solid 
line and denoted as A. The thin lines are labeled with a letter denoting the model, for description see Table Q] We also present the 
sensitivities of future gravitational wave experiments experiments. Bottom panel was created for clarity, because models G, Fl and 
F2 are no different form the standard model. 



7 



I. Kowalska et al.: Gravitational wave background from population III binaries 



10" f 1 1 1 1 1 1 1 r 



a 




f G w [Hz] 

Fig. 5. The gravitational wave background form population III stars. Comparison with the results of iRosadol (1201 lh . Solid line 
corresponds to our standard model A. Dotted and dashed lines are showing background from different classes of compact objects 
calculated bv lRosadol (1201 lh . 
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